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Abstract. Dynamical system models of complex biochemical reaction networks are usu- 
ally high-dimensional, nonlinear, and contain many unknown parameters. In some cases 
the reaction network structure dictates that positive equilibria must be unique for all val- 
ues of the parameters in the model. In other cases multiple equilibria exist if and only if 
special relationships between these parameters are satisfied. We describe methods based 
on homotopy invariance of degree which allow us to determine the number of equilibria 
for complex biochemical reaction networks and how this number depends on parameters 
in the model. 
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1. Introduction 



Dynamical system models of complex biochemical reaction networks are usually high- 
dimensional, nonlinear, and contain many unknown parameters. As was shown recently 
in [CTF06], based on the assumption of mass-action kinetics, graph-theoretical properties 
of some biochemical reaction networks can guarantee the uniqueness of positive equilibrium 
points for any values of the reaction rate parameters in the model. On the other hand, 
relatively simple reaction networks do admit multiple positive equilibria for some values of 
the parameters as shown in [CF05, CF06, CTF06]. 

The aforementioned results do not address the dependence of the number of equilibria on 
the parameter values unless there is a unique equilibrium for every set of parameters. Also 
they do not address the general problem of existence of positive equilibria. Here we describe 
methods using degree theory to analyze general biochemical dynamics (not only mass-action 
kinetics). These methods allow us to determine how the number of positive equilibria for 
a complex biochemical reaction network depends on the parameters of the model. They 
will often also imply the existence of positive equilibria. Also we obtain uniqueness of 
positive equilibria in various situations under significantly weaker assumptions than in 



1.1. Overview. We are interested in equilibria for high- dimensional, nonlinear dynamical 
systems that originate from chemical dynamics. These dynamical systems are systems of 
ordinary differential equations of the form 



where / is a smooth function defined on a subset of the orthant M>o of vectors c in 
having nonnegative components. Such dynamical systems usually have a large number of 
state variables, i.e., n is large. In addition, the parameters defining / are often not well 
known. The focus of this paper is on equilibria for dynamical systems of the form (11. ip . 
that is on c* for which /(c*) = 0. 

We consider the dynamical system (11.11) on a subset Vt of M>o which is the closure of 
a domain VL in M"q. We give conditions for the number of equilibria of (II. ip to remain 
constant as we "continuously deform" (homotopy) the function / through a family of 
functions. A key assumption is that the following condition holds for all members of the 
family: 

(DetSign) The determinant of the Jacobian matrix |f(-) of f is either 
strictly positive or strictly negative on Q. 

(Recall that the Jacobian %{c) at c is the matrix {^(c), i,j = 1, . . . ,n}.) 



[CF05,CF06,CTF06]. 
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What we observe is that when the condition (DetSign) holds for all f in the family and Q is 
bounded, then the number of equilibria for the dynamical system U.l\) is a constant for all f 
in the family, provided there are no equilibria on the boundary of Q for any f in the family 
(see Theorem \l.l\) . We further indicate how this result extends to unbounded domains such 
as Q = ]R"o under suitable conditions (see e.g.. Theorem including those associated 
with a mass- conserving reaction network operating in a chemical reactor with inflows and 
outflows. 



This paper extends previous findings in several ways. 



The (DetSign) condition was introduced by Craciun and Feinberg [CF05,CF06] in the con- 
text of chemistry with Q = R"q and they observed that many chemical reaction networks 
have the property (DetSign). They gave many examples and many tests for this condition 
to hold in the case where / is a system of polynomials and Q is M"q. Then they [CF05,CF06] 
proved that if the components of / are polynomials corresponding to mass-action kinetics 
(operating in a continuous flow stirred tank reactor), and if (DetSign) holds on ]R"q for all 
positive "rate constants", then for each particular choice of rate constant, when an equi- 
librium exists, it is unique. Here we obtain stronger conclusions with weaker assumptions. 
In particular the following are features of our approach. 



(1) Rather than all positive "rate constants" we can select a (vector valued) rate con- 
stant fco of interest at which (DetSign) holds. Then one merely needs a continuous 
curve k{X) of "rate constants" joining ko to another ki at which (DetSign) holds 
and at which the dynamical system has a unique positive equilibrium. 

(2) For a mass-conserving reaction network operating in a chemical reactor with in- 
flows and outflows, under the (DetSign) assumption in (1), we prove existence and 
uniqueness of a positive equilibrium, see Theorem 15. 8[ 

(3) In (1) and (2), the function / need not be polynomial, but is required only to be 
continuously differentiable. Of practical importance are rational / as one finds in 
Michaelis-Menten or Hill type chemical models, see §3, §3 

(4) We give methods, see ^ combining the items above to describe large regions of 
rate constants where a chemical reaction network has a unique positive equilibrium. 



We also point out in this paper that the biochemical reaction network models introduced 
and analyzed by Arcak and Sontag [AS06,AS08] satisfy (DetSign) and we can also rule out 
boundary equilibria (where they give enough data). Consequently, under extremely weak 
hypotheses, we obtain that each of these models has a unique positive equilibrium, see §2.21 
The findings of Arcak and Sontag are impressive in that under strong hypotheses they prove 
global asymptotic stability of equilibria, a topic that this paper does not address. 



1.2. More detail. Now we give some formal definitions. Let f2 be a domain in M", i.e., 
an open, connected set in R". We denote the closure of Q hj n and the boundary of f2 by 
do,. A function / : f2 — M" is smooth if it is once continuously differentiable on Q. If Q 
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is bounded, for such a smooth function /, the following norm is finite: 

||/||o:=sup||/(c)|| 

Here || ■ || denotes the Euclidean norm on M". When Q is bounded, a family fx:Q-^ M" 
for A G [0, 1], is a continuously varying family of functions provided each fx is smooth 
and the mapping A ^ /a is continuous on [0, 1] with the norm || ■ \\^ on the functions fx- 
A zero of / : ^ M" is a value c G such that /(c) = 0, where is the zero vector in 
M". A zero of / is an equilibrium point for the dynamical system (11. ip . 

The following is an immediate consequence of Theorem 13.21 which will be proved in ^ 
This theorem and examples given in ^ are designed to illustrate our approach; then more 
targeted theorems are given in ^and §3 followed by more examples in ^ 

Theorem 1.1. Suppose Q G MJ^ is a bounded domain and /a : ^ M" for A G [0, 1], is 
a continuously varying family of smooth functions such that fx does not have any zeros on 
the boundary of Q for all A G [0, 1]. //det (^(c)) 7^ for all c G Q whenever A = and 
whenever X = 1, then the number of zeros of fo in Q equals the number of zeros of fi in Q. 

The domain of interest for chemical dynamics (cf. (11.11) ) is typically the orthant M>0' 
but this is not bounded, so it violates the hypothesis "fi is bounded". Thus in applying 
Theorem 11.11 we must approximate ]R"o by a large bounded domain Q and check for the 
absence of boundary equilibria. One can think of the boundary dfl in two pieces: that 
which intersects the boundary of M^q, called the sides of Q, and the outer boundary, 
dQ n IR>o- ^® show that if we assume conservation of mass (e.g., by atomic balance) in 
our model and augment with suitable "outflows", then natural bounded domains Q can 
be chosen which have no equilibria on the outer boundary. An example of such a natural 
bounded domain in ]R?,q is shown in Figure [H Also under assumptions of positive invariance 
and augmentation with "inflows", there are no equilibria on the sides. In these cases we 
conclude that there is exactly one nonnegative equilibrium c* for (11.11) and that it is actually 
a positive equilibrium, i.e., it lies in M"q. This result is described in detail in §l]and §3 

I. 3. Organization of the paper. In this paper we give many examples of widely vary- 
ing types to illustrate our contention that our method applies broadly and is easy to use. 
Section [2] gives several examples illustrating Theorem II. 1[ Section [3] summarizes degree 
theory since our proof is based on this and relies on the observation that when (DetSign) 
is true, and there are no boundary equilibria, then "the degree of / with respect to 0" 
equals ± the number of equilibria in a bounded domain. Then in ^ we prove Theorem 

II. 11 and more. Section H] describes the mathematical benefits of mass dissipation (including 
mass conservation) and "inflows" and "outflows" . Section [5] describes a chemical reaction 
network framework which contains in addition to mass-action kinetics, Michaelis-Menten 
and Hill dynamics. We conclude in ^with more examples and new methods presented in 
the context of these examples. For many of our examples, the determinant of the Jacobian 
of / was computed symbolically using Mathematica. As a complement to this paper, we 
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C2 



A side of Q 



T 

A side of 

Figure 1: An example of a natural bounded domain in M?,q with the outer boundary and 
sides indicated. 



have established a webpage at 'http://www.math.ucsd.edu/~helton/chemjac.html' contain- 
ing Mathematica notebooks for many examples in this paper and a demonstration notebook 
that readers may edit to run their own examples. 



2. Examples 

Our goal in this section is to present some examples showing how to use Theorem 1 In the 
process we mention that all chemical reaction examples of Arcak and Sontag [AS06,AS08] 
satisfy (DetSign) and fit well into our approach here. Later in §H]we give broader categories 
of examples. 

2.1. Two examples on treating boundary equilibria. We start with two examples, 
the study of which goes back to a class of examples studied by Thron [T078,T91]. Here, 
c satisfies (11.11) and the Jacobian for all c has the form: 

—ai • • • —bn 
bi —a2 ' ■ • 

62 -0-3 '■• : 

! ■•. ■•. ■•. 

• • • bn-1 -an _ 

where ai > 0, bi > 0, i = 1, . . . ,n, may depend on c. This cyclic feedback structure is 
common in gene regulation networks, cellular signaling pathways, and metabolic pathways 
[AS08]. Thron showed that all eigenvalues of ^ have nonnegative real part (local stability) 




Outer boundary 



(2.1) 



dc 
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if < (sec(7r/n))". Arcak and Sontag showed that an equihbrium of such a dynamical 

system is unique and globally stable under strong global restrictions. 

Here we observe that our key assumption, the determinant of the Jacobian never changes 
sign, is met, which is a major step toward checking when a unique positive equilibrium 
exists for this class of problem. 

Lemma 2.1. When the Jacobian has the form /12. we have 
which if not zero has sign independent of ai,bi > 0. 



(2.2) ^^M^j = (-ir [n^=ia, + n^=,6,, 



Proof: Direct computation. ■ 

We shall now consider several examples from papers of Arcak and Sontag primarily to 
illustrate that checking for absence of boundary equilibria is straightforward; subsequently 
we obtain existence and uniqueness of an equilibrium. In this section we assume that 
all parameters in the reactions are strictly positive. In the sequel, we shall use c as an 
abbreviation for the time derivative of c. 

Example 2.2. We start with an example from §6 of [AS06] which they took from Thron 
[T91]. For this, 

(2.3) Cl 

(2.4) C2 

(2.5) C3 

Now we apply Theorem 11.11 to obtain the conclusion: 



Pi Co 



P3C1 



P2 + C3 
P3C1 - P4C2 

P5C3 

P4C2 - 



n« + Cq 



For each set of parameters pj > 0, j = 1, 2, . . . , 6, Cq > 0, there is a unique 
equilibrium point c* in M^q for the chemical reaction network with dynamics 
given by l^2.3\} - [275\} . and there is no equilibrium point on the boundary of 



However, we have made no comment on stability (even local), while [AS06] gives certain 
conditions ensuring global stability. (In fact, one can algebraically solve for the two equi- 
libria of fi2.3p - fi2.5i) as functions of the parameters. Inspection reveals that exactly one 
of these is in R^q and neither is on the boundary of M>o- The point of us treating this 
example is to show in a simple context how our method works.) 



Proof: We shall apply Theorem 1 1.1 1 to prove that there is a unique equilibrium, inside any 
sufficiently large box; hence there is a unique equilibrium in the orthant. The right hand 
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side fp^coic) of the differential equations fl2.3l) - fl2.5p . while a function of c, also depends on 
positive parameters (p, Cq). One can check that the Jacobian for any of these parameters 
has the form in (12. ip for all c G IR>05 since n = 3 and all parameters are strictly positive, 
the Jacobian determinant is strictly negative. Note that for any two values of the positive 
parameters, (p*, Cq) and (p^ cj), ,,.)+(i_;,)(pt,ct)' ^ ^ 1]' defines a continuously varying 
family of smooth functions on any bounded subset of ffi>o- We check below that the no 
equilibria (i.e., no zeros of /(p,co)) on the boundary hypothesis holds on any sufficiently 
large box, for all positive parameters (p, cq), and thereby conclude using Theorem 1 1 . 1 1 1 hat 
the number of equilibria of fl2.3l) - fl2.5p in R"q does not depend on (p, cq) provided the 
parameters are all strictly positive. Computing the equilibria at one simple "initial" value 
of (p*, Cq) then finishes the proof. 

No equilibria on the boundary of the orthant: Suppose an equilibrium has C2 = 0. Then 
equation fl2.4p implies Ci = which contradicts fl2.3p . Likewise if we start by assuming 
Ci = we get ci > and a contradiction. On the other hand if C3 = 0, then (12. 5p implies 
C2 = 0, which reverts to the case considered first. Thus there are no equilibria on the 
boundary of M>q. 

No equilibria on the outer boundary of some big box: Suppose < 5 < | and ^ > pj > 6 
for all j and Cq < 1. Pick i7 to be a box 1] := {c G M?,o • < {^Y for j = 1,2,3}. An 
equilibrium on the outer boundary of the box satisfies 

(1) ci = (|)^ which by (Q implies {^f > = pgCi > {^f. A contradiction. 
OR 

(2) C2 = (1)"^ which by (!23|) implies | > = P4C2 > (|)^. A contradiction. 
OR 

(3) C3 = (|)^ which by adding (Q, UM, implies that 

^3 ^ j_ > PiCo ^ P5C3 ^ S ■ ji ^ S 

ji ~ P2 + C3 P6 + C3 ~ I + 53 + 1 ■ 

A contradiction. 



Initializing: Up to this point. Theorem 11.11 tells us that each choice of parameters yields the 
same number of equilibria! It is easy to compute for oneself that there is a simple choice of 
parameters which yields a unique positive equilibrium, for example pj = 1 for all j yields 
the unique equilibrium, Ci = C2 = y+^' C3 = Cq. Thus there is one and only one equilibrium 
in M5,q for each value of the positive parameters (p, Cq) . I 

Example 2.3. In Example 1 of §4 in [AS08], the authors describe a simplified model 
of mitogen activated protein kinase (MAPK) cascades with inhibitory feedback, proposed 
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biCi di{l — Cl) fi 



in [K00,SHYDWL01]. For this, 

(2-6) Cl = , ^, , 

Cl + ai Cl + (1 - Cl) 1 + kc3 

fO 7^ ■ ^2C2 d2il - C2) 

2.7 C2 = ■ + — rCi 

C2 + 02 e2 + (1 - C2) 

(2.8) = 'f-''\ '-: 

C3 + 03 Cs + I 1 - C3) 



The variables Cj G [0,1], j = 1,2,3 denote the concentrations of the active forms of 
the proteins, and the terms 1 — Cj, j = 1,2,3, indicate the inactive forms (after non- 
dimensionahzation and assuming that the total concentration of each of the proteins is 1). 
Here the parameters ai, 02, 03, 61, di, ci, 62, (^2, 62, &3, d^, 63, /i, k are strictly positive. 

Let := {c G : < Cj < 1, j = 1, 2, 3} denote the open unit cube, the domain where 
this model holds. Now we show how to apply Theorem 11.11 on Q to conclude that: 



There is a unique equilibrium in Q for any choice of the strictly positive 
parameters, ai, 02, 03, 61, di, Ci, 62, c?2, 62, &3, c?3, 63, /i, fc- 



Proof: First the Jacobian has the form (12.11) . Thus (12.21) implies that the determinant is 
strictly negative for all strictly positive parameters and concentrations c G fi. The proof 
follows the same outline as Example 12.21 Now we check the required items: 

No equilibria on the boundary of the unit cube: Suppose there is an equilibrium c on 
the boundary of Q. Then the equilibrium equations imply that 



(1) If Cl = then (12. 6p forces Ci = 1. Contradiction. 

(2) If Cl = 1 then (12. 6p forces Ci = 0. Contradiction. 

(3) If C2 = then (12. 7p forces Ci = 0. Contradiction as above. 

(4) If C2 = 1 then (12. 7p forces C2 = 0. Contradiction. 

(5) If C3 = then (12. 8p forces C2 = 0. Contradiction as above. 

(6) If C3 = 1 then (12. Sp forces C3 = 0. Contradiction. 



Initializing: [AS08] proves that there are choices of parameters compatible with this model 
for which there is a unique stable equilibrium point in Q. Alternatively, one can compute 
for a simple choice of parameters that there is a unique positive equilibrium. 

The discussion exactly as before implies that there is a unique equilibrium for all fixed 
strictly positive values of the parameters. ■ 
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We mention here that the question of how one finds good intiahzations for the rate constants 
might be a topic for further research. The goal would be to find methods for systemati- 
cally selecting rate constants that produce systems whose equilibria can be determined by 
analytic means. We have not explored this topic at all. 



2.2. The theory of Arcak and Sontag. Now we shall make some general comments 
on [AS06,AS08]. There were four chemical reaction examples presented in the two papers 
[AS06,AS08]. So far we have treated two of the four here in this section. The third example, 
Example 2 in §4 of [AS08], is a small variant of Example 12.31 above and it can be treated 
in a similar manner to that example. In particular, it has a Jacobian of the form (12. ip . We 
now turn to the fourth example of Arcak and Sontag. 

Example 2.4. This is Example 3 in §4 of [AS08] which we do not describe in detail, since 
it requires about a page. While its Jacobian does not have the form (12. ip . it is easy to 
analyze (using Mathematica) and what we found is that the determinant of the Jacobian 
of / is positive at all strictly positive c. Thus the theory described here applies provided 
suitable boundary behavior holds. Boundary behavior was not possible to determine since 
the example was a rather general class whose boundary behavior was not specified. In a 
particular case where more information is specified one might expect that this could be 
done. ■ 

Arcak and Sontag [AS08] present a general theory which contains the examples considered 
in this section and which does not match up simply with ours. Their theory assumes 
an equilibrium exists (we do not). It places global restrictions on the equilibrium which 
guarantee that it is a unique globally stable equilibrium (we address uniqueness but not 
stability). However, while the Arcak-Sontag theory is different than ours, we do point out 
in this section that all four of their chemical examples have Jacobians whose determinant 
sign does not depend on chemical concentration, so our approach applies directly, and with 
a bit of attention to boundary behavior, gives existence of a unique positive equilibrium. 
However, we do not obtain the very impressive global stability in [AS08]. 



3. Degree and homotopy of maps 

The proof of Theorem 1 1 . 1 1 and other results in this paper is based on classical degree theory. 
The degree of a function is invariant if we continuously deform (homotopy) the function 
and we use that to advantage in this paper. 

Now we give the setup. If f2 C M" is a bounded domain, and if a smooth (once continuously 
differentiable) function f : Q has no degenerate zeros, and has no zeros on the 

boundary of Q, then the topological degree with respect to zero of f (or simply the degree 
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of f) equals 




where sgn : M — > { — 1, 0, 1} is the sign function, Zf is the set of zeros of / in Q, and c* is 
a degenerate point means det (c*)) = 0. The de gree of a map naturally extends from 
nondegenerate smooth functions to continuous functions f : Q ^ M". For this, one can 
approximate / uniformly with smooth functions Fk that have no degenerate zeros and no 
zeros on the boundary of Q, and then define the degree of / as the limit of the degrees of 
Fk. The key fact is: this construction of the degree of / is independent of the approximates 
Fk. Fortunately, we shall only need to compute deg(/) on smooth nondegenerate /. For a 
quick account of this theorem and the main properties of degree, see Ch 1.6A of [B77]. 

Homotopy invariance of the degree is the following well known property: 

Theorem 3.1. Consider some bounded domain C M" and a continuously varying family 
of smooth functions /a : ^2 — >■ for A G [0, 1], such that f\ does not have any zeros on the 
boundary ofQ for all X G [0, 1]. Then deg{fx) is constant for all X G [0, 1]. 

Now we give a slightly more general theorem than Theorem 11.11 stated in the introduction. 

Theorem 3.2. Suppose Q and fx, X G [0, 1], are as in Theorem \3.1i Then for any X G [0, 1] 

such that det (^(c)) 7^ for all c G il, the number of zeros of fx in must equal the 
absolute value of the degree of fx in Q, which equals the absolute value of the degree of fx' 
for any X' G [0, 1]. 

Proof: If A G [0, 1] is such that the determinant det{dfx/dc) does not vanish in Q, then 
sgn {det{dfx/dc)) is independent of c. This implies that the zeros of fx are nondegenerate 
and, by the formula for the degree of fx, that | deg(/A)| equals the number of zeros of fx 
in Q. The fact that the degree does not vary with A is immediate from Theorem 13. 1[ ■ 

Remark 3.3. For | deg(/A)| to count the number of zeros of fx in fi, sgn (det (^(c*))) need 
only be the same for all zeros c* in Q, not for all c E fl. Sadly this weakening of hypotheses 
is hard to take advantage of in practice. 

Remark 3.4. From the viewpoint of numerical calculation. Theorem 13.21 strongly suggests 
that if the no boundary zeros hypothesis holds, and (DetSign) holds for / = /a at one 
value of A = Ai, and if one can calculate all zeros of fx at some other value of A = A2 
where (DetSign) also holds, then we can determine the number of zeros at A = Ai. Indeed, 
often we can find a A2 for which fx^ is "simple" in the sense that all zeros for A2 are 
non-degenerate and the zeros can be readily computed along with the Jacobians there, 
and consequently deg(/A2) can be computed. The import for numerical calculation is that 
finding a single equilibrium is often not so onerous. After finding one equilibrium one 
typically makes a new initial guess and tries to find another. Knowing if one has found 
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all of the equilibria is the truly daunting task, since it is nearly impossible to ensure this 
by experiment. Thus theoretical results (hopefully those here) help with this very difficult 
computational question. 



4. Mass-dissipating dynamical systems 



In this section we consider a general dynamical system model which includes the more 
specific dynamics of conservative chemical reaction networks, augmented with infiows and 
outfiows, as described in the next section. In chemical engineering, the latter is commonly 
refered to as dynamics that goes with a continuous fiow stirred tank reactor (CFSTR). 
In biochemistry, one may view this as a model for intracellular behavior with production 
and degradation, or with infiow and outfiow across the cell boundary. Here all species 
components are subject to infiow and outfiow, however, to approximate the conservation 
of some species such as enzymes, one may take the associated infiow rate value in Cin and 
degradation factor in Ao to be arbitrarily small, if desired. 

In preparation for defining a dynamical system on the orthant K>o, we consider a smooth 
function g : M"q — > MJ^, where g has the property that for each j G {l,...,n}, the j^^ 
coordinate of g{c) is nonnegative whenever the j^^ coordinate of c G M>o zero. Consider 
the dynamical system associated with this function given by 



(4.1) c = g{c) for c G 



This dynamical system fl4.ip is called positive-invariant because of the condition on g. 
This guarantees that the dynamics leaves the orthant ]R>o invariant. Given m G IR"q, the 
dynamical system (14. ip is called mass-dissipating with respect to m if 

(4.2) m ■ g{c) < 

for all c G M"q; it is called mass- dissipating if it is mass-dissipating with respect to m for 
some m G IR>o- ^^^^ case, on M>0' 

d(m ■ c) 

4.3 A^=^.^c<0. 

dt 



Now we consider the dynamical system (14. ip augmented with infiows and outfiows: 
(4.4) c = Cin- AoC + g{c), 

where A,, is an n x n diagonal matrix with strictly positive entries on the diagonal. We 
interpret the term Cj„ G M"q as a constant inflow rate, and the term AqC as an outflow 
rate which for each component is proportional to the concentration of that component. It 
is easy to check that with this augmentation, the dynamics still leaves the orthant M>o 
invariant. However, the mass-dissipating property is only inherited at large values of the 
concentration c. 
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We are now ready to state our main theorem in this context. 

Theorem 4.1. Let Cin, m G M"q, and Aq be an n x n diagonal matrix with strictly positive 
diagonal entries. Consider a smooth function g : M>o such that the dynamical system 

( 1^. j[ ) is positive-invariant and mass- dissipating with respect to m. Define 

f{c) := dn- AoC + g{c) /or c G M>o. 

Then the augmented system A4-4\ l> with inflows and outflows, has no equilibria on the 
boundary o/M>q, and if det (|^) ^ on M"q, then there is exactly one equilibrium point 
for this system in ]R"q. 



Proof. It suffices to prove that / has no zeros on the boundary of M% and if det (|f ) =^ 
on M"q, then / has exactly one zero in M"q. 

Define 

/a(c) := Q„ - AoC + Xg{c), for c G M%, A G [0, 1]. 
Fix M > m ■ Cin and let 

Qm = {ceRlo-.m- (Aoc) < M}. 

Then is a bounded domain and {fx : A G [0, 1]} is a continuously varying family of 
smooth functions on CIm- For j = 1, . . . ,n, consider c^ G Qm such that the j^^ coordinate 
of is zero. Then the j^^ coordinate of fx{c^) must be strictly positive, because the j^^ 
coordinate of Cj„ is strictly positive, and the j^^ coordinate of g{c^) is nonnegative, by 
the positive- invariance assumption. Therefore f\ has no zeros on the sides of flM, i-e., on 
^Im n dRyQ. Also, we have 

(4.5) m ■ /a(c) = m ■ Cin — rn ■ (Aoc) + Am ■ g{c) < m ■ Cin — m ■ (Aqc) < 

for all c E Qm such that m- (Aqc) = M. Here we have used the mass-dissipating property of 
m for the ffist inequality and the fact that M > m- Cin for the second inequality. It follows 
that fx has no zeros on the outer boundary of Qm, i-e., on {c G M"o • ' (^qc) = M}. 
Thus, fx has no zeros on the boundary of Qm for all A G [0, 1]. Then, by Theorem 13. ![ 
the degree of fx on Qm, deg{fx,^M), is constant for all A G [0,1]. Next we observe 
that c* = (Ao)~^Cj„ is the unique zero of /o and is inside Qm, and ^ = — Aq, and so 
by fl3.ll) . we obtain deg(/o,fiA/) = sgn(det(— Ao)) = (—1)"- Hence, by Theorem 13.21 if 
det (1^) = det (^) 7^ on Qm, then f = fi has exactly one zero in Qm- 

Since M > m ■ Cin was arbitrary and the sets Qm '■ M > m ■ Cin fill out ffi>o, it follows that 
/ has no zeros on the boundary of R^,. If furthermore, det (f ) ^ on all of Rl,, then it 
follows that / has exactly one zero in M^q. ■ 

Remark 4.2. A special case of mass-dissipating is mass-conserving, namely m-g{c) = 0. For 
dynamical systems associated to chemical reaction networks this has a natural interpreta- 
tion. Indeed, the dynamics of chemical concentrations resulting from chemical interactions 
among several types of molecules will be mass- conserving whenever there exists a mass 
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assignment for each chemical species which is conserved by each reaction, or whenever each 
chemical species (or molecule) is made up of atoms that are also conserved by each reaction. 
More generally, the dynamics will be mass-dissipating whenever no reaction produces more 
mass than it consumes, respectively, produces more atoms than it consumes. Mass conser- 
vation implies det(||) = 0, since m- 1| = when m- g = 0. Thus augmenting with outflows 
is required to make the hypothesis on the sign of det(|^) in our theorems meaningful. The 
paper [HKG08], which builds on the current one, introduces a more general determinant 
that applies when there are no outflows (or only some outflows) . This then helps one count 
equilibria in a manner generalizing what we have done here. 

Remark 4.3. Theorem l4.1l still holds with a much less restrictive definition of mass-dissipating, 
e.g., by replacing "m" with the gradient "VL" for an appropriate class of functions L : 
]R>Q IR>o. Here mass or atom count is behaving like what is called storage function in 
engineering systems theory, see [K02]. Indeed the inequality m ■ fx{c) < m ■ Cin — m ■ (Aqc) 
derived in (14.51) is what is called a dissipation inequality on the storage function c —>■ m ■ c, 
which in fact also plays the role of a "running cost" . 



5. Dynamics of chemical reaction networks 



We now introduce the standard terminology of Chemical Reaction Network Theory (see 
[HJ72, F72, F95, CF05]). A chemical reaction network is usually specified by a finite set of 
reactions ^ involving a finite set of chemical species 

For example, a chemical reaction network with two chemical species Ai and A2 is schemat- 
ically given in the diagram 



(5.1) 2Ai , - Ai + A2 . - ■2A2 

The dynamics of the state of this chemical system is defined in terms of functions CAi{t) 
and CAiit) which represent the concentrations of the species Ai and A2 at time t. The 
occurrence of a chemical reaction causes changes in concentrations; for instance, whenever 
the reaction Ai + A2 ^ 1A\ occurs, the net gain is a molecule of Ai, whereas one molecule 
of A2 is lost. Similarly, the reaction 2A2 —>■ 2Ai results in the creation of two molecules of 
Ai and the loss of two molecules of A2. 

A common assumption is that the rate of change of the concentration of each species 
is governed by mass-action kinetics [HJ72, F72, F79, F87, F95, SOI, CF05, CF06, CTF06], 
i.e., that each reaction takes place at a rate that is proportional to the product of the 
concentrations of the species being consumed in that reaction. For example, under the 
mass-action kinetics assumption, the contribution of the reaction Ai + A2 ^ 2Ai to the 
rate of change of c^i has the form kAi+A2^2AiCAiCA2y where kAi+A2~^2Ai is a positive number 
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called the reaction rate constant. In the same way, the reaction 2A2 2Ai contributes 
the negative value —2k2A2-*2Aic\.^ to the rate of change of c^j. Collecting the contributions 
of all the reactions, we obtain the following dynamical system associated to the chemical 
reaction network depicted in (5.1): 



(5.2) C^i = —k2Ai^Ai+A2CAi + kAi+A2^2AiCAiCA2 — k Ai+ A2-^2A2C A^C A2 

+ k2A2'^Ai+A2(^\2 ~ 2A;2Ai-^2A2C^i + 2A;2A2-+2AiC^2 
CA2 = ^2Ai-+Ai+A2CAi — kAi+A2->2AiCAiCA2 + ^Ai+A2->2A2 C^jCAa 
- k2A2~*Ai+A2C% + 2k2Ai-^2A2c\x " 2A;2A2-*2Ai 



The objects on both sides of the reaction arrows (i.e., 2Ai,Ai + A2, and 2^2) are called 
complexes of the reaction network. Note that the complexes are non-negative integer com- 
binations of the species. On the other hand, we will see later that it is very useful to think 
of the complexes as (column) vectors in M", where n is the number of elements of J^, via an 
identification of the set of species with the standard basis of R", given by a fixed ordering of 

,'2 

the species. For example, via this identification, the complexes above become 2Ai 



Ai + A2 



and 2y42 







We can now formulate a general setup which includes 



many situations, certainly those above. 



5.1. The general setup. Now we present basic definitions and illustrate them. 

Definition 5.1. A chemical reaction network is a triple (^5^,^, i^), where =5^ is a set of n 
chemical species, ^ is a finite set of vectors in R>q with nonnegative integer entries called 
the set of complexes, and C ^ x ^ is a finite set of relations between elements of 
denoted y ^ y' which represents the set of reactions in the network. Moreover, the set 3? 
cannot contain elements of the form y ^ y; for any ?/ G ^ there exists some G ^ such 
that either or ?/' —> y; and the union of the supports of all ?/ G ^ is J^, where the 

support of an element a G M" is supp(a) = {j : aj 7^ 0}. To each reaction y ^ y' & 3?, we 
associate a reaction vector given by y' — y. 



The last two constraints of the definition amount to requiring that each complex appears 
in at least one reaction, and each species appears in at least one complex. For the system 
(15. ip . the set of species is ^ = {Ai, A2}, the set of complexes is ^ = {2Ai, Ai + A2, 2A2} 
and the set of reactions is = {2Ai ^ A1 + A2, A1 + A2 ^ 2^2, 2^2 ^ 2Ai}, and consists 
of 6 reactions, represented as three reversible reactions. 
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In examples we will often refer to a chemical reaction network by specifying only, since 
^ encompasses all of the information about the network. In the sequel we shall sometimes 
simply say reaction network in place of chemical reaction network. 

Definition 5.2. A kinetics for a reaction network (=5^,^,=^) is an assignment to each 
reaction y ^ y' & ^ of a. reaction rate function 

By a kinetic system, which we denote by (^5^, K), we mean a reaction network taken 

together with a kinetics. 



For each concentration c e M"q, the nonnegative number Ky^yi{c) is interpreted as the 
occurrence rate of the reaction y ^ y' when the chemical mixture has concentration c. 
Hereafter, we suppose that reaction rate functions are smooth on M>o, and that Ky^yi{c) — 
whenever supp(?/) ^ supp(c). Although it will not be important in this article, it is natural 
to also require that, for each y ^ y' E M the function Ky^yi is strictly positive precisely 
when supp(|/) C supp(c), i.e., precisely when the concentration c contains at non-zero 
concentrations those species that appear in the reactant complex y. 

Definition 5.3. The species formation rate function for a kinetic system {y^'^^M^K) is 
defined by r : M>o ~^ where 

r(c)= Ky^y'{c){y'-y) for c G M 

The associated dynamical system for the kinetic system [y, K) is 

(5.3) c = r(c)= Ky^y'ic){y'-y), 

where c e M>q is the nonnegative vector of species concentrations. 



The interpretation of r(-) is as follows: if the chemical concentration is c G ffi>o, then rj(c) 
is the production rate of species j due to the occurrence of all chemical reactions. To see 
this, note that 

and that y'^ — yj is the net number of molecules of species j produced with each occurrence 
of reaction y — > y'. Thus, the right hand side of the equation above is the sum of all 
reaction occurrence rates, each weighted by the net gain in molecules of species j with each 
occurrence of the corresponding reaction. Note that rj{c) could be less than zero, in which 
case |rj(c)| represents the overall rate of consumption of species j. 
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5.1.1. Special case: Mass-action kinetics. 

Definition 5.4. A mass-action system is a quadruple {^,'^,J^,k), where is a 

chemical reaction network and k = {ky^yi) is a vector of reaction rate constants, so that 
the reaction rate function Ky_,yi : M"q — > M>o, for each reaction y y' & is given by 
mass-action kinetics: 

n 

Ky^yl{c) = ky_y'C^ W 6 = ]^ Cf . 

1=1 

(Here we adopt the convention that 0° = 1.) The associated mass-action dynamical system 
is 

(5.4) c= Yl ky^y'c^y'-y). 



In the vector equation (15.41) . the total rate of change of the vector of concentrations c is 
computed by summing the contributions of all the reactions in Each reaction y y' 
contributes proportionally to the product of the concentrations of the species in its source y, 
that is, c^, and also proportional to the number of molecules gained or lost in this reaction. 
Finally, the proportionality factor is ky^yi. For example, we can rewrite (15. 2p in the vector 
form (15. 4p as 



(5.5) 



2A2 



2 



-1 
1 



+ kAi+A2~^2AiClC2 



+ ^yli+yl2->2A2ClC2 



k2Ai- 



2 

-2A2C1 



+ k2A., 



2 

*2AiC2 



5.2. Mass conservation. Now we see in terms of this setup how one obtains mass con- 
servation as defined in §11 

Definition 5.5. The stoichiometric suhspace S* C M" of a reaction network (^, ^) is 
the linear subspace of spanned by the reaction vectors y' —y , for all reactions y ^ y' & ^. 



Note that, according to (15.31) . for a given value of c, the vector c is a linear combination of 
the reaction vectors. This implies that each stoichiometric compatibility class (co + S')nM"Q 
is an invariant set for the dynamical system (15. 3p with initial condition cq G M>q. 

Definition 5.6. A reaction network (^,^, ^) is called conservative if there exists some 
positive vector m G M"q which is orthogonal to all its reaction vectors, i.e., 

m-iy' -y) =Q 

for all reactions in Then m is called a conserved mass vector. 

Remark 5.7. Each trajectory of a conservative reaction network is bounded. A conservative 
reaction network can have no inflows or outflows (see the definition of inflow and outflow 
in the next section). 
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5.3. Main results. We now consider conservative reaction networks augmented with in- 
flow and outflow reactions (for each of the species). An inflow reaction is a reaction of the 
form — > A and an outflow reaction is one of the form A — > 0, where A is a species. The 
reaction vector y' — y associated with an inflow reaction for species j is the vector containing 
all zeros, except that it has a one in the j^^ position. The reaction vector associated with 
an outflow reaction for species j is the negative of the vector associated with an inflow 
reaction for that species. Here, for the kinetics associated with the inflows and outflows, 
we assume that the reaction rate function for each inflow reaction is a positive constant 
and the value of the reaction rate function for each outflow reaction is a positive constant 
times the concentration of the species flowing out. The latter corresponds to degradation 
of each species at a rate proportional to its concentration. The following theorem may be 
used to determine the number of equilibria for conservative reaction networks augmented 
by such inflows and outflows. It requires a positive determinant condition and is the analog 
of Theorem 4.1 in this context. 

Theorem 5.8. Consider some conservative reaction network {S^,'io,M) with conserved 
mass vector m. Let K be a kinetics for this network with associated species formation rate 
function g. Consider an augmented kinetic system {y, J^, J^) obtained by adding inflow 
and outflow reactions for all species so that the associated dynamical system is: 

(5.6) c = r{c) := Cin - AoC + g{c), 

where q„ G M"o Ao is an n x n diagonal matrix with strictly positive diagonal entries. 
Suppose that 



for all c G ffi>o- Then the dynamical system ( f5. 6\) has exactly one equilibrium c* in ]R"q 
and no equilibria on the boundary o/M>q. 

Proof. We want to apply Theorem 14.11 The function g is given by the right member of 
(15. 3p . where the functions Ky^yi are all smooth and have the property that Ky^y/{c) = 
whenever supp(?/) ^ supp(c). Consequently, g is smooth and, whenever c G M>o is such 
that Cj = for some j, then we have 



because Ky^yi{c) = whenever yj > and Cj = 0, by the support property of K. It 
follows that the dynamical system (15. 6p is positive-invariant. Furthermore, the system is 
mass-dissipating, since 



by the assumption that the reaction network {S^, M) is conservative. The conclusion 
then follows immediately from Theorem 14.11 ■ 




gj{c)> Ky^yic){-yj) =0, 
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Remark 5.9. The results described above use the assumption that all species have inflows, 
in order to conclude that there are no equilibria on the boundary of IR"o- other hand, 

for very large classes of chemical reaction networks described in [ADS07], this assumption 
is actually not needed in order to rule out the existence of such boundary equilibria (an 
observation by David Anderson University of Wisconsin [A]). 

Remark 5.10. The paper [BDB07], for the case of "nonautocatalytic reactions", gives a 
condition on the "stoichiometric matrix" (in our terminology the matrix whose columns 
are the vectors y'—y for y y' & which is necessary and sufficient for the determinant of 
the Jacobian of r to be of one sign for all concentrations and all Ky^yi which are monotone 
increasing in each variable. Our theory is less restrictive as illustrated by Example 16.11 

6. Applications 

In practice, most dynamical system models of biochemical reaction networks contain a large 
number of unknown parameters. These parameters correspond to reaction rates and other 
chemical properties of the reacting species. In this section we treat a variety of examples of 
such models and illustrate the use of Theorems 14. II and 15. 8[ In some of these examples, the 
determinant of the Jacobian det (|^) is of one sign everywhere on the open orthant for all 
parameters and in some it is not. Even in the latter cases we describe ways to find classes 
of rate functions for which there exists a unique positive equilibrium. 

The first subsection assumes mass-action kinetics and defines (reminds) the reader of the 
Craciun-Feinberg "determinant expansion" via an example. Critical is the sign of each 
term in the expansion and whether all terms have the same sign or miss this by "a little" , 
namely, only one or two terms in the determinant expansion has a coefficient with an 
anomalous sign. Here we point out that all examples in [CF05,CF06,CTF06] have at most 
one or two anomalous signs. When there are no anomalous signs, these papers show that 
any positive equilibrium is unique for all parameter values, and they develop and use graph 
theoretical methods for determining when there are no anomalous signs. Here, for cases of 
few anomalous signs, we propose and illustrate a technique for identifying parameter values 
for which a positive equilibrium exists and is unique. The paper, [HKG08], subsequent to 
this one, gives ways of counting the number of anomalous signs in terms of graphs associated 
to a chemical reaction network. 

In the second subsection, we continue with the general framework of §5l and move be- 
yond mass-action kinetics to rate functions satisfying certain monotonicity conditions (see 
Definitions 16.31 and 16.41) . The weaker condition. Definition 16. 4[ holds for many biochem- 
ical reactions and allows one to make sense of the signs which occur in the determinant 
expansion. Hence one can apply the methods in this paper. 

The number of anomalous signs can be determined for the examples in this section using: 
(a) the graph-theoretic methods of Craciun and Feinberg [CF05, CF06] when there are no 
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anomalous signs and the kinetics are of mass-action type, and (b) symbolic computation of 
the determinant of the Jacobian using Mathematica when there are some anomalous signs 
or the kinetics are general (not necessarily mass-action). The reader will find Mathematica 
notebooks at 

http: / /www. math. ucsd.edu/~helton/ chemjac.html 

for many of the examples in this section (including all of those that fall under (b)), as 
well as a demonstration notebook that readers may edit to run their own examples. This 
software works well when the number of species is small; for larger numbers, the determinant 
expansion has too many terms to be handled readily. 

We conclude this preamble with some intuition underlying the case when there are anoma- 
lous signs. In general, we expect that for very small values of the parameters appearing in 
the reaction rate functions for a conservative reaction network (and, in the limit, for vanish- 
ing parameter values), the dynamics of the system augmented by inflows and outflows will 
be dominated by the inflow and outflow terms, and (\.ei{dr / dc) will not vanish; moreover, if 
the inflow and (linear) outflow terms dominate the dynamics, then the equilibrium will be 
unique. Examples 16. II and 16.61 illustrate how this observation can be made rigorous and can 
be used together with Theorem 15.81 and the proof of Theorem 14. II to conclude the existence 
and uniqueness of an equilibrium for a subset of the parameter space, even if the result 
does not hold for the entire parameter space. 



6.1. The determinant expansion, its signs and uses. 

Example 6.1. Consider the mass-action kinetics system given by the chemical reaction 
network (16.1 1) , which is an irreversible version of the network shown in Table 1.1 (i) of [CF05] 
(see Table lU^i) below): 

A + B P 

(6.1) B + C ^ Q 

C ^ 2A 

If we add inflow and outflow reactions for all species, the associated dynamical system 
model for (16.11) is 



Ca 


= fco- 


-^A 


- kA- 


^qCA — kA+B- 


^pCaCb + 2kc- 


-*2ACc 


Cb 


= ^0- 


^B 


~kB- 


^qCb — kA+B- 


^pCaCb - kB^ 


-C^qCbCc 


Cc 


= ^0- 


■^C 


-kc- 


■*oCc - kB+c- 


-^qCbCc - kc- 


*2ACc 


Cp 


= h- 


■*P 


— kp- 


^OCp + kA+B- 


^PCaCb 






= ^0- 




- kq- 


^qCq + kB+C- 


^QCbCc- 





According to Remark 4.3 in [CF05] the dynamical system above does have multiple positive 
equilibria for some values of the reaction rate parameters. 
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If we assume that all outflow rate constants /ca^O; ■■■} ^q^o are equal to 1, then the deter- 
minant of the Jacobian of the reaction rate function is: 

det{dr/dc) = -1 - kA+B^pCA - ks+c^QCc - kB+c^QCs 

— kB+C^qkA+B^pCACB — kc^2A — kc^2AkA+B^pCA 
(6.3) - kc^2AkB+C^QCc - kA+B^pCB - kA+B-*pkc^2ACB 

— kA+B^pkB+C^Q^B " kA+B^pkB+C^QCBCc 
+ kA+B-.pkB+C-^Qkc^2ACBCc- 

Note that there is only one positive monomial in the expansion in (16. 3p . The concentrations 
in it are cbCc, but there is also a negative monomial with concentrations cbCc, and the 
two combine to give 

[- kA+B~,pkB+C~*Q + kA+B^pkB+C^Qkc^2A(^BCc- 

Thus if kc^2A < 1, then the positive monomial will be dominated by a negative monomial. 
Therefore, if kc^2A < l, then det (9r/9c) ^ for this network, everywhere on ]R5,g. 

Note that (m^, m^, mc, mp, mq) = (1, 1, 2, 2, 3) is a conserved mass vector for the reaction 
network (16. ip . It follows from Theorem 15.81 that (16. 2p . the dynamical system for (16. ip . 
augmented with inflows and outflows (with outflow rate constants equal to one), has a 
unique positive equilibrium for all positive values of the reaction rates such that kc^2A < 1. 
Note that this uniqueness conclusion would not follow directly from the theory of [CF05, 
CF06] nor from that in [BDB07], since these works pertain only when the determinant has 
the same sign for all rate constants and species concentrations. 

The same method can be applied to conclude that the reversible version of the reaction 
network (16.11) . augmented with inflows and outflows (with outflow rate constants set equal 
to one), also has a unique positive equilibrium for all positive values of the reaction rates 
such that kc^2A < 1; moreover, even if the (positive) outflow rate constants fc^^o, kq^o 
are not necessarily equal to 1, the same conclusion holds if we know that kc^2A < kc^o- 
■ 

Example 6.2. Here we summarize several examples with mass-action kinetics (in the next 
subsection we consider some of these reactions with more general kinetics). Of the eight 
examples in [CF05,CF06], which are chemical reaction networks augmented with inflows 
and outflows (with outflow rate constants equal to one), two have the property that the 
coefficients of the terms in their Jacobian determinant expansion all have the same sign, 
and the other six have all but one sign the same. The first observation is from [CF05] and 
the second observation, emphasizing that there is only one anomalous sign, is new here. An 
analysis as in Example 16.11 can be applied here. Table [1] is a list of the examples showing 
how many "anomalous" signs each determinant expansion has. 

A similar accounting holds for examples of reaction networks in [CTF06], see Table 1, 
page 8699. These reactions involve enzymes which [CTF06] treat with mass-action kinetic 
models. They find reaction networks 1,2,3,5,7,9 in this table do not have any anomalous 
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Reaction 

llcL W(Ji Jv 


Num. of "anomalous" signed 

LcllllD 111 LlCt t;-A.|JcLilt)lUil 


(i) A + B^P B + C 


1 


(ii) A + B^P B + C^Q 





C + D ^ R D + E ^ S 

E:^2A 


1 


(iv) A + B^P B + C 
C^A 





(v) A + B^F A + C^G 
C+D^B C+E^D 


1 


(vi) A + B ^ 2A 


1 


(vii) 2A + B^'dA 


1 


(viii) A + 25 ^ 3A 


1 



Table 1: Some examples of reaction networks and the signs of coefficients in their Jaco- 
bian determinant expansion when augmented with inflows and outflows (with outflow rate 
constants equal to one). 



signs. Here we point out that the remaining reaction networks, 4, 6 and 8 have very few 
anomalous signs. The reaction network 4 is 

S + E^ES-^E + P, I + E^EI, I + ES ^ ESI ^ EI + S 

and has only 1 "anomalous" sign, and the reaction network 6 is 

Sl + E^ESl, S2 + E^ES2, S2 + ES1 ^ ES1S2 ^ S1+ES2, ES1S2^E + P 

and has only 2 "anomalous" signs. Reaction network 8 has 4 anomalous signs out of a 
total of over 3000 terms. Here all reactions are augmented by outflows with outflow rate 
constants set to one. (For the cases of no anomalous signs these outflow rate constants can 
be taken to be arbitrarily small without changing the answeiQ.) The theory of Sections [Hand 
|5] applies, if there are (arbitrarily small) inflows and outflows, to yield that there is a unique 
positive equilibrium, for reaction networks 1,2,3,5,7,9. It also leaves open the possibility 
that one can apply the technique in Example 16.11 to get a unique positive equilibrium for 
certain rate constants in reaction networks 4, 6, 8. These applications of our theory require 
finding a conserved "mass" for the system without infiow and outfiow, which is easy to do 
in all cases. 

6.2. General reaction rate functions. In this subsection, we follow the setup in ^and 
move beyond mass-action kinetics to a very general classes of rate functions. 



See [CF06iee] for how one can eliminate outflows for the enzymes. 
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Definition 6.3. We say that a reaction rate function Ky^yi is consumptively increasing, 
if for each species i belonging to the support of the vector y, the partial derivative of the 
reaction rate function, dKy^y/ /dci, is strictly positive on the open orthant. 

It is very common to assume that the reaction rate functions Ky^yi are consumptively 
increasing, since this simply means that the rate of a reaction increases whenever the con- 
centration of a consumed species is increased unilaterally. In particular, the consumptively 
increasing property is true for many common chemical reaction rate laws, such as many 
Michaelis-Menten and Hill laws, as well as for all mass-action kinetics [KS98]. All examples 
in this section entail consumptively increasing reaction rates. 

The consumptively increasing property can fail to hold for some classes of reactions in- 
cluding those involving inhibitory enzymes and for those in which a Michaelis-Menten rate 
depends on the products of the reaction [Rec81]. However, the next more lenient assump- 
tion handles these and many additional biochemical situations. 

Definition 6.4. We say that a reaction rate function Ky^y' is strictly monotone, if for 
each species i on which the function Ky^yi actually depends, the partial derivative of the 
reaction rate function, dKy^yi / dci, has one strict sign on the open orthant. 



More generally, the main technique used in this section is to compute the determinant 
expansion of the Jacobian as a sum of terms, each of which is a product of partial derivatives 
dKy^yi / dci where i belongs to the support of y. For strictly monotone rate functions we 
can assign a ± to each term according to whether that term is everywhere positive or 
negative on the domain M"q. That is, strict monotonicity guarantees the technique of 
tracking anomalous signs in the determinant expansion applies. 

Examples 12.21 12.31 and 12.41 which involve inhibitory feedback can be written in the form 
(15. 3p with strictly monotone rate functions. As we observed in ^the determinant of the 
Jacobian, (|^), is positive for all of these situations. However, at this point the terminology 
is in place so that we can mention the more refined property that each of these examples 
has no anomalous signs. 

Example 6.5. Consider the chemical reaction network (16. 4p . which is the reversible network 
shown in (ii) in Tabled] but, unlike in [CF05], in this example we don't assume that the 
kinetics is mass-action. 

A + B ^ P 

(6.4) B^C ^ Q 

C + D ^ R 
D ^ 2A 

We augment this reaction with inflows and outflows where the outflow matrix Ao is nor- 
malized to be the identity. We suppose that each of the reaction rate functions Ky_,yi is 
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consumptively increasing as in Definition 16.31 We compute the expansion of det{dr / do) in 
terms of the partial derivatives dKy^y'{c)/dci, for i belonging to the support of y. It is 
a sum of coefficients times monomials in these partial derivatives; the set of coefficients is 
shown in fl6.5l) . 

{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 



(6.5) 



-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-2, 


-2, 


-2, 


-2, 


-2, -2, 


-2, 


-2, 


-2, 


-2, 


-2, 


-2, 


-2, 


-3, 


-1, 


-1, 


-1, 


-1, 


-1,-1, 


-1, 


-3, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1,-1, 


-1, 


-2, 


-2, 


-2, 


-2, 


-2, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1,-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-2, 


-2, 


-2, 


-2,-2, 


-2, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-2, 


-2, 


-1, 


-1,-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-2, 


-2, 


-2, 


-2, 


-2,-2, 


-2, 


-2, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-2, 


-2, 


-2,-1, 


-1, 


-1, 


-1, 


-1, 


-1, 


-2, 


-2, 


-2, 


-1, 


-1, 


-1, 


-2, 


-1} 





To summarize this list there are 138 terms in the expansion of det{dr/dc) and the set 
of coefficients of these terms contains exactly: 96 minus ones, 40 minus twos, 2 minus 
threes and no positive terms. This is more information than we need, since the point is 
that all these numbers are negative. This implies that det{dr/dc) ^ for this network, 
for all values of c G M^g "consumptively increasing" rate laws. Also, note that 

{mA,mB, mcrriD, mp, niQ, mji) = (1, 1, 1, 2, 2, 2, 3) is a conserved mass vector for the reac- 
tion network (16. 4p . Therefore, we can apply Theorem 15.81 to conclude that for the reaction 
network (16.41) with rate laws which are all consumptively increasing, and with inflows and 
outflows (with rate constants equal to one), has a unique positive equilibrium. ■ 

Example 6.6. This example is exactly parallel to Example 16.51 except that now we consider 
the chemical reaction network which is the network shown in (v) in Table [TJ Unlike in 
[CF05] , in this example we don't assume that the kinetics is mass-action. We assume that 
the reaction rate functions are consumptively increasing. Also we augment with inflows 
and outflows where the outflow matrix A^ is normalized to be the identity matrix. We find 
that there are 167 terms in the expansion of det{dr/dc) involving the partial derivatives 
of dKy^yi I dci for i belonging to the support of y and the set of coefficients of these terms 
contains exactly: 146 minus ones, 20 minus twos and one positive term. The positive term 
is 

K'b^C+d{cb) K'o^c+e{cd) ^i+C-.G(cA, Cc) i^A+B^F(cA, Cb) 

Here F*-^'"^ (resp. F'^'^'^^) denotes the partial derivative of F with respect to its first (re- 
spectively second) variable. 

One can think of many conditions on the reaction rate functions that make the determinant 
have one sign on the open orthant. Typically, the more complicated they look, the more 
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lenient is the assumption. Here are some examples. We first note that the reaction network 
(v) of Table [1] is mass-conserving with mass vector m = (1, 3, 1, 2, 1, 4, 2). 

(1) We can collect all terms containing K'j-f^(j_^_j^[cD) K^^q^q{ca, cc) K^^^_^p[cA, cb), 
this yields (-1 + K's^c+DicB)) K'j, 

Thus if we assume that 1 > K'^^(^^j^{cb) for all cb > 0, then the determinant is 
negative on the open orthant. 

(2) Alternatively, we can collect all terms containing K^A+c^ci'^A, cc) K^A+B^pi'^A, cb) 
and extract its coefficient which yields 

(6.6) - [1 - K'^^cM^b)] K'^^cMod) - [1 + K^ciLBicc, cd)] [1 + 4+Lij(cc, cn)]. 

Thus assuming this is negative for all positive cb,cc,ce) makes the determinant 
negative on the open orthant. We see that the second requirement is less stringent 
than the first. 

(3) For any M > m ■ the boundary of the set Qm '■= {c G M^q '■ ^ - c < M} contains 
no equilibria, as in the proof of Theorem 14. 1[ A yet weaker assumption than (1) 
and (2) is that the function (16.61) is negative on Qm for a particular M > m ■ Cm- 

In cases (1) and (2), we can apply Theorem 15.81 to conclude that for the reaction network 
in Table [T] (v) with rate laws which are all consumptively increasing, after augmentation 
with infiows and outfiows (where Aq is the identity matrix), there exists a unique positive 
equilibrium. For case (3) we can apply the proof of Theorem 14.11 to conclude under the 
same conditions that there is a unique equilibrium in fijv/ . B 

Remark 6.7. We emphasize that the homotopy-based methods described in this paper not 
only imply uniqueness, but also existence of a positive equilibrium for many dynamical 
systems derived from chemical reaction networks, while the methods in [CF05,CF06] only 
imply uniqueness of an equilibrium. Also, as we saw, our methods may be used for models 
containing very general chemical kinetics laws (not necessarily mass-action kinetics). 
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